import numpy as np

import matplotlib.pyplot as plt

import time

import math

returnlist = list()

#Fills array with value specified, with potential for half split.
def fill(ls, val, split):

	if split != val:
		for a in range(len(ls)/2):
			ls[a] = val
		for b in range(len(ls)/2, len(ls)):
			ls[b] = split
	else:
		for a in range(len(ls)):
			ls[a] = val
	
	return ls

#Returns a list of specified <size> with value 0.0.
def emptygrid(size):

	ls = list()

	for a in range(size):
		ls.append(0.0)

	return ls

#displacement has three portions:
#Past = [0]
#Present = [1]
#Future = [2]
#The past and present are either initial value determined, or
#are derived from the array of the next time step. The future timestep
#is determined from past and present steps, the result of both steps becomes a
#rolling record of the past two steps.

def mainloop(time, xspacing, tspacing, modulus=1.0, othermodulus=1.0, howdense=1.0, drivingfreq=1.5):

	#Setup
	size = int(round(10/xspacing))

	displacement = list()

	density = howdense
	modulus = fill(emptygrid(size), modulus, othermodulus)
	
	
	x = [x*xspacing for x in range(size)]
	
	upast = emptygrid(size)

	upresent = emptygrid(size)

	ufuture = emptygrid(size)

	displacement.append(upast)
	displacement.append(upresent)
	displacement.append(ufuture)
	
	changelist = emptygrid(size)

	courantcond(max(modulus), density, xspacing, tspacing)

	#Interate

	for a in range(int(round(time/tspacing))):

		step(modulus, density, displacement, a, xspacing, tspacing, size, drivingfreq)	

	storage(displacement)
	
	#Plot

	plt.axis([0, 10, -2, 2])

	plt.plot(x, displacement[1])

	plt.show()


def step(mod, density, displacement, curtime, xspacing, tspacing, size, drivingfreq):

	shake(displacement[1], curtime, tspacing, drivingfreq)
	
	ischanging = checkmod(mod)

	#d^2(u)/dx^2 = d/dt(du/dt)
	for pos in range(1, size - 1):
		if True:
			displacement[2][pos] = changingmod(mod, density, pos, displacement[:2], xspacing, tspacing)
		else:	
			displacement[2][pos] = paulsolve(mod, density, pos, displacement[:2], xspacing, tspacing)
		
	displacement[0] = displacement[1][:]
		
	displacement[1] = displacement[2][:]
	

def courantcond(modulus, density, xspacing, tspacing):
	result = abs(math.sqrt(modulus/density))*tspacing/xspacing

	if result > 1:
		print('Warning, these values exceed the courant condition!!')
		print('Your result could be unstable!')	

#Determines if 
def checkmod(mod):

	diffval = mod[0]
	for a in mod:
		if a != diffval:
			return True
		else:
			diffval = a		
	return False	

def storage(curlist):

	returnlist = curlist[:]

	alist = curlist[:] 


	def showlist(val):

		print(alist[val])


#Simulates a forced motion at the end of the rod
def shake(u, curtime, tspacing, drivingfreq):

	sin, pi = np.sin, np.pi
	
	u[0] = sin(drivingfreq*curtime*tspacing)

#Approximates d^2(u)/dx^2 (depricated)
def paulsolve(mod, density, pos, displacement, xspacing, tspacing):

	val = (mod[pos]/density)*(tspacing**2.0/xspacing**2.0)*(-2.0*displacement[1][pos] + displacement[1][pos-1] + displacement[1][pos+1]) + 2.0*displacement[1][pos] - displacement[0][pos]

	return val

#Function calculates numerically calculated time step in displacement.
def changingmod(mod, density, pos, displacement, xspacing, tspacing):

	du = ((displacement[1][pos]-displacement[1][pos-1])/xspacing+(displacement[1][pos+1]-displacement[1][pos])/xspacing)/2.0

	dE = ((mod[pos]-mod[pos-1])/xspacing+(mod[pos+1]-mod[pos])/xspacing)/2.0
	
	val = 2.0*displacement[1][pos] - displacement[0][pos] + (tspacing**2.0/density)*dE*du + (mod[pos]/density)*(tspacing**2.0/xspacing**2.0)*(-2.0*displacement[1][pos]+displacement[1][pos-1]+displacement[1][pos+1])

	return val	
